Self-Organized Criticality and Universality in a Nonconservative Earthquake Model 
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We make an extensive numerical study of a two dimensional nonconservative model proposed 
by Olami-Feder-Christensen to describe earthquake behavior. By analyzing the distribution of 
earthquake sizes using a multiscaling method, we find evidence that the model is critical, with no 
characteristic length scale other than the system size, in agreement with previous results. However, in 
contrast to previous claims, we find convergence to universal behaviour as the system size increases, 
over a range of values of the dissipation parameter, a. We also find that both "free" and "open" 
boundary conditions tend to the same result. Our analysis indicates that, as L increases, the 
behaviour slowly converges toward a power law distribution of earthquake sizes -P(s) ~ with 
exponent r ~ 1.8. The universal value of r we find numerically agrees quantitatively with the 
empirical value (r = B + 1) associated with the Gutenberg-Richter law. 
PACS numbers: 05.40.+j, 91.30.Px 



I. INTRODUCTION 

Earthquakes may be the most dramatic example of 
self-organized criticality (SOC) that can be observed 
by humans on earth. Most of the time the crust of the 
earth is at rest, or quiescent. These periods of stasis 
are punctuated by sudden, thus far unpredictable bursts, 
or earthquakes. According to the empirical Gutenberg- 
Richter (GR) law the distribution of earthquake 
events is scale-free over many orders of magnitude in en- 
ergy. The GR scaling extends from the smallest measur- 
able earthquakes, which are equivalent to a truck passing 
by, to the most disastrous that have been recorded where 
hundreds of thousands of people have perished. 

The relevance of SOC to earthquakes was first pointed 
out by Bak and Tang Q, Sornette and Sornette as 
well as Ito and Matsusaki According to this the- 
ory, plate tectonics provides energy input at a slow time 
scale into a spatially extended, dissipative system that 
can exhibit breakdown events via a chain reaction pro- 
cess of propagating instabilities in space and time. The 
GR law arises from the system of driven plates building 
up to a critical state with avalanches of all sizes. These 
above-mentioned authors used a spatially extended (but 
conservative) cellular automata model as a prototype re- 
sembling earthquake dynamics which gave a power law 
distribution of avalanches, or earthquakes. This was fol- 
lowed by a study using block-spring models by Carl- 
son and Langer ||^, who found characteristic earthquake 
sizes, rather than asymptotic critical behavior. Later 
studies of a continuous "train" block-spring model by de 
Sousa Vieira ^ recovered criticality. The train model 
describes a driven elastic chain sliding over a surface 
with friction. It has been conjectured to be in the same 
universality class as interface depinning and a model of 
avalanches in granular piles, which agrees with numerical 
simulation results jl^ . 

Several groups made lattice representations of the 



block-spring model. These models were nonconservative 
|]lT| , ^ and were driven uniformly, but did not display 
SOC. Then Olami, Feder, and Christensen (OFC) intro- 
duced a nonconservative model on a lattice that displayed 
SOC In this simplified earthquake model, sites on 

a lattice are continuously loaded with a force. After a 
threshold force is reached, the sites transfer part of their 
force to their local neighborhood when discharging. Each 
discharge event is accompanied by a local loss in accumu- 
lated force from the system, when the force on each ele- 
ment is reset to zero. A uniform driving force is slowly ap- 
plied to all the elements and the model is completely de- 
terministic. This conceptually simple and seemingly nu- 
merically tractable model reproduces some of the qualita- 
tive phenomenology of the statistics of earthquake events 
such as power law behavior over a range of sizes, inter- 
mittency or clustering of large events [|l3|, and lack of 
predictabihty 0. 

Although SOC and this type of modelling approach has 
been more or less accepted as a reasonable description of 
the phenomena of earthquakes (see for example Ref. jl^] 
and references therein) , the OFC model itself has had a 
controversial existence, both on the theoretical [l^-[l9|] 
and numerical side |r^ , |20| - |2^ , p^ , p5| . The initial numeri- 
cal studies found that the distribution of earthquake sizes 
obeyed finite size scaling (FSS) over the range of system 
sizes that could be studied at the time [^Sj. This placed 
the nonconservative model into the framework of stan- 
dard critical behavior. However, these simulations also 
indicated that there was no universality. In particular 
the exponents characterizing the power law distributions 
appeared to vary with both the dissipation parameter, 
a, and the form of boundary condition. If this were the 
case then the OFC model would be very different from 
familiar critical systems where most microscopic details 
are irrelevant and have no effect on critical coefficients. 
In fact an argument was made [ p^ that one should not 
expect universal behavior in far from equilibrium criti- 
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cal phenomena. If this were correct, it would drastically 
limit the application of any known theoretical tools to 
these problems. 

Another strange aspect was that the dimension D char- 
acterizing the scaling of the cutoff in the earthquake size 
distribution was found numerically to be larger than two. 
This is inconsistent with the fact that each site can only 
discharge a finite number of times in an earthquake event, 
requiring D < 2 for the two dimensional model . This 
last result together with the strange lack of universality 
cast some doubt on whether the OFC model was actually 
critical or just close to being critical, with some large as 
yet undetermined length scale beyond which the earth- 
quake distribution would always be cut off. Hwa and 
Kardar as well as Grinstein and collaborators postulated 
that conservation of the quantity being transported was 
required for criticality |26| , but the theoretical arguments 
made do not take into account SOC phenomena such as 
avalanches and long-term memory associated with the 
self-organization process (for more details see |^^). The 
fact that the random neighbor version of the nonconser- 
vative OFC model is never critical but has an essential 
singularity as the conservative limit is approached [ p^Jl9| 
has added to the mystery surrounding the behavior of the 
lattice model. 

In a previous large scale numerical simulation study of 
the model discussed here, Grassberger Q also claimed 
that the model was critical but found that some of the 
conclusions of OFC "have to be modified considerably." 
He argued that the probability distribution of earthquake 
sizes does not show ordinary FSS over the larger range 
of systems he was able to study, and "conjecture(d) that 
the cutoff of P{s) becomes a step function for L oo," 
although he did not present direct numerical evidence of 
this. 

There are a number of important, unresolved questions 
about the behavior of the model, which have enormous 
implication for any type of eventual theoretical under- 
standing. Is the nonconservative model on a lattice (or 
for fixed connectivity matrix) critical? If so, is the criti- 
cal behavior of the model universal over a range of values 
of a, or for different boundary conditions? Is it described 
by a power-law distribution at all? Are there any other 
universal quantities? What type of data analysis tech- 
nique besides FSS would be useful to extract the large 
scale behavior of the nonconservative model? Our nu- 
merical study and analysis will address those issues and 
answer those questions. 



A. Summary 

In the first section we review the definition of the OFC 
model and present some numerical data for the distri- 
bution of earthquake sizes using standard FSS. For the 
range of lattice sizes we have simulated (up to linear size 



L = 512), our data confirm the lack of apparent FSS 
in the model, particularly in the cutoff region. In Sec- 
tion III, we present an extensive set of results using a 
multiscaling method. We analyze the rescaled probabil- 
ity distribution, p3| , in terms of the quantity 
Dav = logs/ log L, with s being the size of an earth- 
quake. We observe that there are no avalanches with 
Dav larger than two, consistent with the bound imposed 
on the cutoff Sco (see previous discussion). By analyzing 
how this distribution behaves for different values of the 
nonconservation parameter, a, and system size, L, we 
show how the multiscaled probability distribution tends 
to converge to a universal curve as L increases. The direc- 
tion of convergence on increasing L changes as a varies 
enabling us to put fairly firm limits on the asymptotic 
curve. The model appears not be to described at all by 
FSS. However, for s < Sco the distribution converges to- 
ward a power law with a universal exponent r ~ 1.8 over 
a range of a values. Moreover the cutoff in this distribu- 
tion becomes very sharp as L increases and its behavior 
indicates that Sco const(a)L^ as L — > oo. In Section 
IV we summarize our main conclusions. 



II. DEFINITION OF THE MODEL 

We consider a two-dimensional square lattice of L x L 
sites. At each site i, a force Fi is assigned to be a real 
variable. Initially, the force at each site is chosen ran- 
domly from the uniform distribution between and 1. 
The dynamics proceeds by two steps in the limit of infi- 
nite time scale separation between the slow drive, repre- 
senting motion of the tectonic plates, and the earthquake 



process |13 



1. Increase the force at all sites: Find the largest force 
Fmax in the system and increase the force at all 
sites by the same amount 1 — Fmax- 

2. Relax all unstable sites, i.e. sites with Fi > 1: The 
force of an unstable site is reset to zero, Fi 0, 
and a fraction of it, a-Fi, is distributed to each of 
its four nearest neighbours, F„„ Fn„-\-aFi. This 
step is repeated in a parallel update until there are 
no unstable sites left. 

This two step rule is iterated indefinitely. The sequence 
of toppling events (step 2) between application of the uni- 
form drive (step 1) defines an avalanche or earthquake. 
Since only a fraction, 4a, of the force is redistributed in 
each toppling, the model is nonconservative for a < 1/4. 

To completely define the model we need to specify the 
boundary conditions, by defining the dissipation param- 
eter, a, at the sites on the boundaries (af,c). As in OFC, 
we consider both "free" and "open" boundary conditions. 
The sites on the boundaries of the system can be con- 
sidered to be bounded by fictitious sites with Fi = — oo, 
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which can never discharge and only absorb force from the 
boundary sites. In the case of open boundary conditions, 
the sites at the boundary have the same a parameter as 
all other sites in the bulk {abc = a). In the case of free 
boundary conditions, the sites at the boundary have their 
a parameters modified in order to have the same level of 
dissipation per reset event as for sites in the bulk. This 
latter condition implies ck{,c = ol/{1 — q), except at corner 
sites where abc.c = ~ 2a). The model with periodic 
boundary conditions is not critical |^,^,^, and will 
not be discussed. It is probably worthwhile to underline 
at this point that the model is completely deterministic, 
the only possible source of randomness coming from the 
initial conditions. 

After a transient period of many earthquakes, the 
model settles into a statistically stationary state. One 
way to characterize this state is to measure statistical 
properties of the earthquakes. The size of an earthquake, 
s, is defined as the number of resets of the local force 
i'i — > in the system in between applications of the uni- 
form force. One can also measure the temporal duration, 
t, in terms of the parallel update, or the radius of gyra- 
tion r of the sites which participated in the earthquake 
event. 



A. Finite Size Scaling 



We focus on the probability distribution of earthquake 
sizes, s, in a system of size L, Pl{s). If the model is criti- 
cal then this distribution will have no scale other than the 
physical extent L and the lattice constant, which is set 
to unity. One ansatz that can describe critical behavior 
is the FSS ansatz that was previously used by OFC: 

P^{s)^L-PG{^) (1) 



where G is a suitable scaling function, and (3 and D are 
critical exponents describing the scaling of the distribu- 
tion function. As shown in Fig. 1, the model does not 
exhibit FSS. In this figure we chose D ~2 as the largest 
possible allowed value. Still the cutoff in the "collapsed" 
probability distribution moves to the right as L increases. 
Nevertheless, for earthquake sizes smaller than the cut- 
off, this figure shows what appears to be a convergence 
to a well-defined power law, Pl{s) ''^ as L increases 
with the power law exponent r = P/D 2± 1.8 for both 
a = 0.18 and a = 0.21, and possibly also for a = 0.15. 
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FIG. 1. Finite-size scaling plots of Pl{s) in systems with 
open boundary conditions for (a) a = 0.15, (b) a = 0.18, and 
(c) a = 0.21. The critical exponents are D — 2 and /3 = 3.6, 
the slope of the straight line is t = 1.8. Statistics are derived 
from at least 10^ avalanches per data set. For visual clarity, 
curves (b) and (c) have been shifted along the horizontal axis, 
X —i X + 1 and x ^ x + 2 respectively. 



III. MULTISCALING ANALYSIS 

The FSS ansatz is only one possible description of crit- 
ical behavior. As pointed out some time ago by Kadanoff 
and co-workers, some SOC phenomena are better de- 
scribed by a multifractal ansatz, rather than FSS [ p8[ . 
This form has recently been used to clarify the behav- 
ior of the Bak-Tang-Weisenfeld ^ model by Stella and 
co-workers, who have measured different moments asso- 
ciated with the distribution |2^. For the OFC model, it 
appears to us that a clearer picture can be obtained by 
simply examining the probability distribution directly. 

The multiscaling ansatz postulates for the probability 
distribution function Pl(s) a form 

log PUs/ So) ^ / logs/5o \ 
log(L//<,) yiogL/lJ ' ^ > 

where So and lo are parameters which typically (but not 
always) reflect phenomena at small scales associated with 
the lattice Usually, a multiscaling analysis consists 
of choosing these two parameters in order to get the best 
collapse for different system sizes. This is quite different 
than FSS where the critical exponents themselves, re- 
flecting behavior at large scales, must be chosen in an ad 
hoc manner in order to obtain the "best" collapse. We do 
not attempt to collapse the data using the multiscaling 
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form of Eq. 2. Instead, wc define Zq = So = 1 to rep- 
resent the smallest earthquake which occurs at only one 
site and involves only one discharge. Moreover, we define 
the dimension of an earthquake of size s in a systems of 
size L as 

Da„ = log s/ log L , 

and we denote with D^ax the largest value of Day 

In Fig. 2 we show the multiscaled probability distri- 
bution according to Eq. 2, for different system sizes for 
a = 0.21. One observes immediately that all avalanches 
have dimension Dav < 2, as required, with the largest 
dimension Dmax approaching two as the system size in- 
creases. Also it is clear that the shape of the cutoff func- 
tion is sharpening as the system size increases. Since the 
largest dimension D„iax cannot be larger than two, we 
can infer from this that the cutoff region narrows to the 
region Dav 2 and becomes increasingly sharp. 
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ing asymptotic term, which we propose is a power law as 
suggested by Fig. 2. This is 



log Pl{s) 
log(L) 



= F{Dav) = -{TDay) + F^utoffiDav) , (3) 



where Fcutof f in the limit L —>■ oo should be constant up 
to a cutoff near Dmax- 
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FIG. 3. Plots of F^^toff for (a) a = 0.15, (b) a = 0.18, 
(c) a = 0.21, and (d) a — 0.23; we set t = 1.8. Boundary 
conditions are open. Different curves correspond, from left to 
right, to L = 32, 64, 128, 256, 512. 



FIG. 2. Multiscaling plot of Pl{s) for a = 0.21 and open 
boundary conditions. The slope of the straight line is r = 1.8. 

Note that the increase of Dmax as L increases is to- 
tally inconsistent with the notion that the OFC model 
is slightly off criticality, because in that case one would 
expect that the relative size of the largest dissipating 
events with respect to the maximum total force allowed 
in the system would decrease in larger systems. In fact 
what happens is exactly the opposite. In larger systems 
a larger fraction of the total energy can be dissipated in 
the largest event that occurs, and Dmax increases with 
L. This result is completely consistent with the non- 
conservative model being critical, rather than slightly off 
criticality. 

In order to get more explicit visual information on the 
probability distribution we try to subtract out the lead- 



We get a consistent picture for a range of a values by 
choosing t = 1.8, as shown in Fig. 3. Although it ap- 
pears for small system sizes that smaller values of a give 
a steeper distribution with a larger value of r (so the 
lines tend to decrease from left to right rather than re- 
maining horizontal), it is clear from this figure that as 
the system size increases all the different values of a ap- 
pear to reach the same value of r = 1.8, corresponding 
to a completely horizontal line in this figure. The devia- 
tion for small systems is more pronounced for a = 0.15, 
and less again for a = 0.18, becoming the smallest for 
a = 0.21. For a closer to the conservative limit, as shown 
for a = 0.23, the approach to the asymptotic horizontal 
behavior changes direction. Namely smaller systems ap- 
pear to have a smaller exponent r than larger systems, at 
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least for small avalanches. Thus rather than having the 
slope increase as L increases, for a >~ .21, the apparent 
slope decreases as L increases, as clearly demonstrated 
in this figure. 

We ascribe the change of direction to a crossover effect 

of the conservative fixed point. For a close to 1/4, smaller 
avalanches behave as avalanches would in the conserva- 
tive system. It is only the larger avalanches that are af- 
fected by nonconservative dissipation. This is associated 
with the fact that as a approaches 1/4, each site can top- 
ple more and more times in a single earthquake. For any 
finite value of dissipation (1 —4,0}, the maximum number 
of times that a site can reset in an earthquake is deter- 
mined by this dissipation and is finite in the limit of large 
L. However for small systems, the largest avalanches are 
not large enough to be effected by dissipation, and the 
cutoff in the number of times that a site can topple is 
not determined by (1 — Aa) but by L. Then the cutoff in 
the number of topplings at a given site is the same as in 
a conservative system of the same size. 

The results we have described thus far are for the 
model with open boundary conditions. Fig. 4 shows that 
the same behavior occurs for the OFC model with free 
boundary conditions, with the same value r = 1.8. In 
this case, the asymptotic behavior as L increases is ap- 
proached for decreasing apparent slope for both a = 0.17 
and a = 0.20, as is plainly evident in the figure. Again 
the cutoff appears to sharpen as L increases and approach 
D — 2 
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FIG. 4. Plots of Fcutoff in systems with free boundary 
conditions for (a) a = 0.17 and (b) a = 0.20. The expo- 
nent has been set to r = 1.8 as before. System sizes are 
L = 32, 64, 128,256,512. 

Our analysis indicates that the OFC model exhibits 
SOC with a universal power law distribution of events 
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for logs < {Dmax^ogL). As shown in 



Fig. 3, the cutoff gets sharper and sharper as L increases 
with Dmax approaching 2 from below as L increases for 
all values of a we have studied. As indicated above, the 
incorrect results obtained by OFC were due to the fact 
that there is a strong system size dependence that varies 
with a and is not described by FSS. In addition to a sys- 
tematic change in apparent t as L increases, the largest 
dimension of earthquakes, Dmax, is increasing only slowly 
towards two, and probability distribution of avalanche 
dimensions is getting sharper at the cutoff. This means 
that large avalanches are suppressed in small systems rel- 
ative to the total amount of force in the system as com- 
pared to a larger system, and that a FSS "fit" for any L 
range will always give an apparent D > 2, which is not 
allowed. In this sense the model appears to violate FSS 
for all values of L. 



IV. CONCLUSIONS 

The main results of this paper are as follows: The 
nonconservative model on a two dimensional lattice self- 
organizes into a critical state. The critical state is robust 
and universal over a range of values of the dissipation pa- 
rameter, a, and for different boundary conditions. The 
model does not exhibit finite-size scaling. The cutoff be- 
comes sharper as L increases with largest earthquakes of 
dimension Dmax approaching two from below. Neverthe- 
less, the probability distribution of earthquake sizes is a 
power law with a universal exponent r c± 1.8. This value 
can be identified with an exponent B = r — 1 for the dis- 
tribution of energy dissipated in earthquakes. According 
to the Gutenberg-Richter law this is a power law with 
B = 0.8 to 1, completely consistent with our result. 

We thank K. Christensen and H. J. Jensen for help- 
ful conversations, and P. Bak for comments on the 
manuscript. S.L. was supported by the EPSRC through 
a post-doctoral research fellowship. 
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